(*

*****************************************
* WRITTEN BY SIJA [ONLY DIFFGEO] - 2024 *
*****************************************

********************************************************************************
* SOME EXAMPLES OF PROCEDURES DEFINITIONS FROM CLASSICAL DIFFERENTIAL GEOMETRY *
********************************************************************************


*********
* Tools *
*********


01. Resize := proc(P,m::posint,n::posint)
02. ExpGraphAC := proc(P, fName, fExt, widthVal, heightVal)
03. ExpGraphRR := proc(P, fName, fExt, widthVal, heightVal)
  

*********************
* Parametric Curves *
*********************


01. Tan_vec3D := proc(X, unon)
02. Nor_vec3D := proc(X, unon)
03. Bin_vec3D := proc(X, unon)

04. curvature_3D := proc(X)
05. curvatureR_3D := proc(X)

06. torshion_3D := proc(X)
07. torshionR_3D := proc(X)

08. curveTanPln_3D := proc(X, u, w)
09. curveTanPlnVec_3D := proc(X, tanVec, norVec, u, w)

10. curveNorPln_3D := proc(X, u, w)
11. curveNorPlnVec_3D := proc(X, binVec, norVec, u, w)

12. curveRecPln_3D := proc(X, u, w)
13. curveRecPlnVec_3D := proc(X, binVec, tanVec, u, w)

14. sphereRadius := proc(X)
15. sphereCenter := proc(X)
16. osculSphere := proc(X)

17. osculCircleRadius := proc(X)
18. osculCircleCenter := proc(X)

19. Circle_3D := proc(circleCenter, circleRadius, v1, v2)

20. rollCircCentInTanPln_3D := proc(X, r, inORout)
21. rollPntOnCircInTanPln_3D := proc(X, r, inORout)

18. rollCircCentInRecPln_3D := proc(X, r, inORout)
19. rollPntOnCircInRecPln_3D := proc(X, r, inORout)



***********************
* Parametric Surfaces *
***********************


20. surXuTanVec_3D := proc(X, varL, unon)
21. surXvTanVec_3D := proc(X, varL, unon)

22. surTanVec_3D := proc(X, varL, unon)
23. surNorVec_3D := proc(X, varL, unon)

24. coeFirFunForm_3D := proc(X, varL)
25. wFirFunForm_3D := proc(X, varL)
 
26. coeSecFunForm_3D := proc(X, varL)
27. coeThiFunForm_3D := proc(X, varL)

28. surTanPln_3D := proc(X, varL, w, s)
29. surNorPlnXu_3D := proc(X, varL, w, s)
30. surNorPlnXv_3D := proc(X, varL, w, s)

31. surTanLine_3D := proc(X, varL, PointL)
32. surNorLine_3D := proc(X, varL, PointL)

33. Christoffel_3D:= proc(X, varL, indexL)

34. surGerGeoEqu_3D := proc(X, varL)

*)



#***********************
#  Tools - Procedures  *
#***********************


Resize := proc(P,m::posint,n::posint)
      op(0,P)(op(P),ROOT(BOUNDS_X(0),BOUNDS_Y(0),
                   BOUNDS_WIDTH(n),BOUNDS_HEIGHT(m)));
end proc:


ExpGraphAC := proc(P, fName, fExt, widthVal, heightVal)
  local ResizeMe, pl3dStrR, ftype, path, currDir, fileName;
  
  ResizeMe := proc(P,m::posint,n::posint)
      op(0,P)(op(P),ROOT(BOUNDS_X(0),BOUNDS_Y(0),
                   BOUNDS_WIDTH(n),BOUNDS_HEIGHT(m)));
  end proc:
  
  pl3dStrR := ResizeMe(P, widthVal, heightVal):
  
  ftype := convert(fExt, symbol);
  
  currDir := interface(worksheetdir):
  fileName := cat(fName,".", fExt);  
  path := FileTools:-JoinPath([currDir, fileName]):
  
  if FileTools[Exists](path)
  then FileTools[Remove](path); end if;
  
  plotsetup(ftype,plotoutput=path);
  print(pl3dStrR);
  pl3dStrR;
  plotsetup(inline);
  
  plotsetup(default);
  fclose(path);
    
end proc:


ExpGraphRR := proc(P, fName, fExt, widthVal, heightVal)
  local ResizeMe, pl3dStrR, ftype, path, currDir, fileName;
  
  ResizeMe := proc(P,m::posint,n::posint)
      op(0,P)(op(P),ROOT(BOUNDS_X(0),BOUNDS_Y(0),
                   BOUNDS_WIDTH(n),BOUNDS_HEIGHT(m)));
  end proc:
  
  pl3dStrR := ResizeMe(P, widthVal, heightVal):
  
  ftype := convert(fExt, symbol);
  
  currDir := interface(worksheetdir):
  fileName := cat(fName,".", fExt);  
  path := FileTools:-JoinPath([currDir, fileName]):
  
  if FileTools[Exists](path)
  then FileTools[Remove](path); end if;

  plotsetup(ftype, plotoutput = path);
  print(pl3dStrR);
  pl3dStrR;
  
  plotsetup(default);
  fclose(path);
    
end proc:


#*********************************
# Parametric Curves - Procedures *
#*********************************


Tan_vec3D := proc(X, var, unon) # X <- Column Vector
  local XL, Xd, XLx_d, XLy_d, XLz_d; 
      
  XL := convert(X, list);
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  # Column Vector of Tan_vec3D
  if unon
  then 
    simplify(eval(Xd/norm(Xd, 2))) assuming var::real;
  else   
    simplify(eval(Xd)) assuming var::real;
  end if;       
end proc:


Nor_vec3D := proc(X, var, unon) # X <- Column Vector
  local XTV, XL, Xdd, XLx_dd, XLy_dd, XLz_dd, vcp, vcpf; 
  
  XTV := Tan_vec3D(X, var, unon);
      
  XL := convert(X, list);
      
  XLx_dd := diff(XL[1], var, var);
  XLy_dd := diff(XL[2], var, var);
  XLz_dd := diff(XL[3], var, var);
  
  Xdd := < XLx_dd, XLy_dd, XLz_dd >;
  
  vcp := VectorCalculus:-CrossProduct(XTV, Xdd);
  vcpf := VectorCalculus:-CrossProduct(vcp, XTV);
  
  # Column Vector of Nor_vec3D
  if unon
  then
    simplify(eval(vcpf/norm(vcpf, 2))) assuming var::real;
  else   
    convert(simplify(eval(vcpf)), Vector) assuming var::real;
  end if; 
        
end proc:


Bin_vec3D := proc(X, var, unon) # X <- Column Vector
  local XTV, XNV, XL, Xdd, XLx_dd, XLy_dd, XLz_dd, vcp, vcpf; 
  
  XTV := Tan_vec3D(X, var, unon);
  XNV := Nor_vec3D(X, var, unon);
   
  vcp := VectorCalculus:-CrossProduct(XTV, XNV);
    
  # Column Vector of Bin_vec3D
  if unon
  then   
    simplify(eval(vcp/norm(vcp, 2))) assuming var::real;
  else
    convert(simplify(eval(vcp)), Vector) assuming var::real;
  end if;       
end proc:


curvature_3D := proc(X, var) # X <- Column Vector
  local XTV, XL, XLx_dd, XLy_dd, XLz_dd, Xdd, vcp; 
       
  XTV := Tan_vec3D(X, var, false);   
  
  XL := convert(X, list);
  
  XLx_dd := diff(XL[1], var, var);
  XLy_dd := diff(XL[2], var, var);
  XLz_dd := diff(XL[3], var, var);
  
  Xdd := < XLx_dd, XLy_dd, XLz_dd >;
  
  vcp := VectorCalculus:-CrossProduct(XTV, Xdd);
  
  # Curvature in 3D of vector function X
  simplify(eval(norm(vcp, 2)/norm(XTV, 2)^3)) assuming var::real;
     
end proc:


curvatureR_3D := proc(X, var) # X <- Column Vector
  local cur3d; 
  
  cur3d := curvature_3D(X, var);
    
  # Curvature Radius in 3D of vector function X
  if cur3d <>0 
  then 
    abs(simplify(eval(1/cur3d))) assuming var::real; 
  else 
    print("The radius of curvature is infinitely large");
    infinity; 
  end if;
     
end proc:


torshion_3D := proc(X, var) # X <- Column Vector
  local XL, Xd, XLx_d, XLy_d, XLz_d,
                Xdd, XLx_dd, XLy_dd, XLz_dd,
                Xddd, XLx_ddd, XLy_ddd, XLz_ddd,
                vcp, dop, vcpd; 
      
  XL := convert(X, list);
   
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  XLx_dd := diff(XLx_d, var);
  XLy_dd := diff(XLy_d, var);
  XLz_dd := diff(XLz_d, var);
  
  Xdd := < XLx_dd, XLy_dd, XLz_dd >;
  
  XLx_ddd := diff(XLx_dd, var);
  XLy_ddd := diff(XLy_dd, var);
  XLz_ddd := diff(XLz_dd, var);
  
  Xddd := < XLx_ddd, XLy_ddd, XLz_ddd >;
  
  vcp := VectorCalculus:-CrossProduct(Xdd, Xddd);
  dop := VectorCalculus:-DotProduct(Xd, vcp);
  vcpd := VectorCalculus:-CrossProduct(Xd, Xdd);
  
  # Torshion in 3D of vector function X
  simplify(eval(dop/norm(vcpd, 2)^2)) assuming var::real;
  
end proc:
 

torshionR_3D := proc(X, var) # X <- Column Vector
  local tor3d; 
  
  tor3d := torshion_3D(X, var); 
    
  # Torshion Radius in 3D of vector function X
  if tor3d <>0 
  then 
    abs(simplify(eval(1/tor3d))) assuming var::real; 
  else 
    print("The radius of torsion is infinitely large");
    infinity; 
  end if;
  
end proc:


curveTanPln_3D := proc(X, var, u, w) # X <- Column Vector
  local tanvec, norvec;
     
  tanvec:= Tan_vec3D(X, var, true);
  norvec:= Nor_vec3D(X, var, true);
      
  # Parametric tangent plane to X 
  simplify(eval([X[1] + u*tanvec[1] + w*norvec[1], 
                 X[2] + u*tanvec[2] + w*norvec[2], 
                 X[3] + u*tanvec[3] + w*norvec[3]])) assuming var::real;
                
end proc:


curveTanPlnVec_3D := proc(X, var, tanVec, norVec, u, w) # X, tanVec, norVec <- Column Vectors functions of the same one variable
     
  # Parametric tangent plane to X 
  simplify(eval([X[1] + u*tanVec[1] + w*norVec[1], 
                 X[2] + u*tanVec[2] + w*norVec[2], 
                 X[3] + u*tanVec[3] + w*norVec[3]])) assuming var::real;
                
end proc:


curveNorPln_3D := proc(X, var, u, w) # X <- Column Vector
  local binvec, norvec;
      
  binvec:= Bin_vec3D(X, var, true);
  norvec:= Nor_vec3D(X, var, true);
      
  # Parametric normal plane to X 
  simplify(eval([X[1] + u*binvec[1] + w*norvec[1], 
                 X[2] + u*binvec[2] + w*norvec[2], 
                 X[3] + u*binvec[3] + w*norvec[3]])) assuming var::real;
                
end proc:


curveNorPlnVec_3D := proc(X, var, binVec, norVec, u, w) # X, tanVec, norVec <- Column Vectors functions of the same one variable
    
  # Parametric normal plane to X 
  simplify(eval([X[1] + u*binVec[1] + w*norVec[1], 
                 X[2] + u*binVec[2] + w*norVec[2], 
                 X[3] + u*binVec[3] + w*norVec[3]])) assuming var::real;
                
end proc:


curveRecPln_3D := proc(X, var, u, w)  # X <- Column Vector
  local binvec, tanvec;
     
  binvec:= Bin_vec3D(X, var, true);
  tanvec:= Tan_vec3D(X, var, true);
      
  # Parametric rectifying plane to X 
  simplify(eval([X[1] + u*binvec[1] + w*tanvec[1], 
                 X[2] + u*binvec[2] + w*tanvec[2], 
                 X[3] + u*binvec[3] + w*tanvec[3]])) assuming var::real;
                
end proc:


curveRecPlnVec_3D := proc(X, var, binVec, tanVec, u, w) # X, tanVec, norVec <- Column Vectors functions of the same one variable
  local binvec, norvec;
  
  # Parametric rectifying plane to X 
  simplify(eval([X[1] + u*binVec[1] + w*tanVec[1], 
                 X[2] + u*binVec[2] + w*tanVec[2], 
                 X[3] + u*binVec[3] + w*tanVec[3]])) assuming var::real;
                
end proc:


sphereRadius := proc(X, var) # X <- Column Vector
  local XL, XLx_d, XLy_d, XLz_d, Xd, cur3D, cur3D_d, 
        tor3d, c_b; 
    
  XL := convert(X, list);
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
           
  cur3D := curvature_3D(X, var);
  cur3D_d := diff(cur3D, var);
  tor3d := torshion_3D(X, var);
  
  simplify(eval(sqrt(1/cur3D^2 + cur3D_d^2/(norm(Xd, 2)^2*cur3D^4*tor3d^2)))) assuming var::real;
 
end proc:


sphereCenter := proc(X, var) # X <- Column Vector
  local XL, XLx_d, XLy_d, XLz_d, Xd, 
        cur3D, cur3D_d, tor3d, XBV, XNV, c_b;
        
  XL := convert(X, list);
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  cur3D := curvature_3D(X, var);
  cur3D_d := diff(cur3D, var);
  
  tor3d := torshion_3D(X, var);
  
  XBV := Bin_vec3D(X, var, true);
  XNV := Nor_vec3D(X, var, true);
  
  # Osculation sphere center
  simplify(eval(X + (1/cur3D)*XNV - cur3D_d/(norm(Xd, 2)*cur3D^2*tor3d)*XBV)) assuming var::real;
 
end proc:


osculSphere := proc(X, var) # X <- Column Vector
  local sC, sR, sp;
  global theta, phi;

  sC := sphereCenter(X, var); 
  sR := sphereRadius(X, var);
  # Attention: theta and phi are globa variables
  sp := convert([cos(theta)*sin(phi)*sR, sin(theta)*sin(phi)*sR, cos(phi)*sR], Vector)   assuming var::real;
    
  # Osculation sphere           
  simplify(eval(sC + sp));
  
end proc:


osculCircleRadius := proc(X, var) # X <- Column Vector
  local cur3d; 
  
  cur3d := curvature_3D(X, var);
    
  # Osculation Circle Radius in 3D of vector function X
  if cur3d <>0 
  then 
    abs(simplify(eval(1/cur3d))) assuming var::real; 
  else 
    print("The radius of the osculation circle is infinitely large");
    infinity; 
  end if;

end proc:


osculCircleCenter := proc(X, var) # X <- Column Vector
  local cR, unvec;
  
  cR:= osculCircleRadius(X, var);
  unvec:= Nor_vec3D(X, var, true); 
 
  # Osculation circle center          
  simplify(eval(X + cR*unvec)) assuming var::real;
     
end proc:


Circle_3D := proc(circleCenter, circleRadius, v1, v2, var)
  local xc, yc, zc, center, radius;
  
  center := circleCenter; 
  radius := circleRadius;
    
  # Parametrize circle
  # s := 0..2*Pi is a circle parametr:
  xc := center[1] + radius*cos(s)*v1[1] + radius*sin(s)*v2[1];
  yc := center[2] + radius*cos(s)*v1[2] + radius*sin(s)*v2[2];
  zc := center[3] + radius*cos(s)*v1[3] + radius*sin(s)*v2[3];
  
  # Cirrcle on a plane spread over vectors v1 and v2 in 3D
  simplify(eval([xc, yc, zc])) assuming var::real;
  
end proc:


rollCircCentInTanPln_3D := proc(X, var, r, inORout) # X <- Column Vector
  local unvec;
  
  unvec:= Nor_vec3D(X, var, true);
   
  if inORout 
  then
    # A circle with this center rolls along the "concavity" - in of the curve
    simplify(eval(X + r*unvec)) assuming var::real;
  else
    # A circle with this center rolls on the "convexity" - out of the curve   
    simplify(eval(X - r*unvec)) assuming var::real;
  end if;
  
end proc:


rollPntOnCircInTanPln_3D := proc(X, var, r, inORout) # X <- Column Vector
  local XL, XLx_d, XLy_d, XLz_d, Xd, tanvec, norvec, rollCC, 
        arcelem, rarc, rotang;
  
  XL := convert(X, list);
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  arcelem := sqrt(VectorCalculus:-DotProduct(Xd, Xd));
  arcelem := subs(var = u, arcelem);
  
  rotang:= evalf(Int(arcelem, u = 0..var)/r) assuming var::real;
      
  tanvec := Tan_vec3D(X, var, true);
  norvec := Nor_vec3D(X, var, true);
  rollCC := rollCircCentInTanPln_3D(X, var, r, inORout);
              
  # Rolling Point On Circle In Tangent Plain in 3d
  if inORout
  then   
    rollCC - r*tanvec*cos(rotang) + r*norvec*sin(rotang);
  else
    rollCC + r*tanvec*sin(rotang) - r*norvec*cos(rotang);
  end if;

end proc:


rollCircCentInRecPln_3D := proc(X, var, r, inORout) # X <- Column Vector
  local ubvec;
  
  ubvec:= Bin_vec3D(X, var, true);;
  
  if inORout 
  then
    # A circle with this center rolls along the "concavity" - in of the curve
    simplify(eval(X + r*ubvec)) assuming var::real;
  else
    # A circle with this center rolls on the "convexity" - out of the curve   
    simplify(eval(X - r*ubvec)) assuming var::real;
  end if;
  
end proc:


rollPntOnCircInRecPln_3D := proc(X, var, r, inORout) # X <- Column Vector
  local XL, XLx_d, XLy_d, XLz_d, Xd, tanvec, binvec, rollCC, 
        arcelem, rarc, rotang;
  
  XL := convert(X, list);
    
  XLx_d := diff(XL[1], var);
  XLy_d := diff(XL[2], var);
  XLz_d := diff(XL[3], var);
  
  Xd := < XLx_d, XLy_d, XLz_d >;
  
  arcelem := sqrt(VectorCalculus:-DotProduct(Xd, Xd));
  arcelem := subs(var = u, arcelem);
  
  rotang:= evalf(Int(arcelem, u = 0..var)/r) assuming var::real;
    
  tanvec := Tan_vec3D(X, var, true);
  binvec := Bin_vec3D(X, var, true);
  rollCC := rollCircCentInRecPln_3D(X, var, r, inORout);
            
  # Rolling Point On Circle In Tangent Plain in 3d
  if inORout
  then   
    rollCC - r*tanvec*cos(rotang) + r*binvec*sin(rotang);
  else
    rollCC + r*tanvec*sin(rotang) - r*binvec*cos(rotang);
  end if;
  
end proc:


#***********************************
# Parametric Surfaces - Procedures *
#***********************************


surXuTanVec_3D := proc(X, varL, unon) # X <- Column Vector
  local XL, varU, XuL, Xu, norXu; 
  
  XL := convert(X, list);             

  varU :=  varL[1];
  
  XuL := diff(XL, varU)  assuming varU::real;
  Xu := simplify(convert(XuL, Vector), assume=real);
  norXu := simplify(eval(Xu/norm(Xu, 2)), assume=real);  
  
  if unon
  then
    norXu;
  else   
    Xu;
  end if;

end proc:


surXvTanVec_3D := proc(X, varL, unon) # X <- Column Vector
  local XL, varV, XvL, Xv, norXv; 
  
  XL := convert(X, list);
  varV :=  varL[2];
  
  XvL := diff(XL, varV) assuming varV::real;
  Xv := simplify(convert(XvL, Vector), assume=real);
  norXv := simplify(eval(Xv/norm(Xv, 2)), assume=real);
  
  if unon
  then
    norXv;
  else   
    Xv;
  end if;

end proc:


surTanVec_3D := proc(X, varL, unon) # X <- Column Vector
  local TanV, varU, varV, norTanV;
  (*
  varU := varL[1];
  varV := varL[2];
  
  assume = [varU::real, varU::real];
  *)
  TanV := simplify(surXuTanVec_3D(X, varL, unon) + surXvTanVec_3D(X, varL, unon), assume=real);
  norTanV := simplify(eval(TanV/norm(TanV, 2)), assume=real);
  
  if unon
  then 
    norTanV;
  else   
    TanV;
  end if;

end proc:


surNorVec_3D := proc(X, varL, unon) # X <- Column Vector
  local Xu, Xv, NorV, varU, varV, norNorV;
  (*
  varU := op(1, indets(X,name));
  varV := op(2, indets(X,name));
  
  assume = [varU::real, varU::real];
  *)
  Xu := surXuTanVec_3D(X, varL, unon);
  Xv := surXvTanVec_3D(X, varL, unon);
  
  NorV := simplify(VectorCalculus:-CrossProduct(Xu, Xv), assume=real);
  norNorV := simplify(eval(NorV/norm(NorV, 2)), assume=real);
  
  if unon
  then 
    norNorV;
  else   
    NorV;
  end if;
  
end proc:


coeFirFunForm_3D := proc(X, varL) # X <- Column Vector
  local Xu, Xv, E, F, G;
  
  Xu := surXuTanVec_3D(X, varL, false);
  Xv := surXvTanVec_3D(X, varL, false);
  
  E := simplify(eval(VectorCalculus:-DotProduct(Xu, Xu)), assume=real);
  F := simplify(eval(VectorCalculus:-DotProduct(Xu, Xv)), assume=real);
  G := simplify(eval(VectorCalculus:-DotProduct(Xv, Xv)), assume=real);
  
  [E, F, G];
  
end proc:


wFirFunForm_3D := proc(X, varL) # X <- Column Vector
  local cFFF, E, F, G;
  
  cFFF := coeFirFunForm_3D(X, varL);
  E := cFFF[1];
  F := cFFF[2];
  G := cFFF[3];
  
  simplify(E*G - F^2, assume=real);
  
end proc:


coeSecFunForm_3D := proc(X, varL) # X <- Column Vector
  local Xu, Xv, L, M, N, varU, varV, unNorV,
        XuL, XvL, XuuL, XuvL, XvvL, Xuu, Xuv, Xvv;
  
  varU := varL[1];
  varV := varL[2];
  
  Xu := surXuTanVec_3D(X, varL, false);
  Xv := surXvTanVec_3D(X, varL, false);
  
  XuL := convert(Xu, list);             
  XvL := convert(Xv, list);
  
  XuuL := diff(XuL, varU)  assuming varU::real;
  XuvL := diff(XuL, varV)  assuming varV::real;
  XvvL := diff(XvL, varV)  assuming varV::real;
  
  Xuu := convert(XuuL, Vector);
  Xuv := convert(XuvL, Vector);
  Xvv := convert(XvvL, Vector);
  unNorV := surNorVec_3D(X, varL, true);
  
  L := simplify(eval(VectorCalculus:-DotProduct(unNorV, Xuu)), assume=real);
  M := simplify(eval(VectorCalculus:-DotProduct(unNorV, Xuv)), assume=real);
  N := simplify(eval(VectorCalculus:-DotProduct(unNorV, Xvv)), assume=real);
  
  [L, M, N];
  
end proc:


coeThiFunForm_3D := proc(X, varL) # X <- Column Vector
  local A, B, C, varU, varV, unNorV, unNorVL,
        unNorVLu, unNorVLv, unNorVu, unNorVv;
    
  varU := varL[1];
  varV := varL[2];
  
  unNorV := surNorVec_3D(X, varL, true);
  unNorVL := convert(unNorV, list);
  
  unNorVLu := diff(unNorVL, varU)  assuming varU::real;
  unNorVLv := diff(unNorVL, varV)  assuming varU::real;
  
  unNorVu := convert(unNorVLu, Vector);
  unNorVv := convert(unNorVLv, Vector);
  
  A := simplify(eval(VectorCalculus:-DotProduct(unNorVu, unNorVu)), assume=real);
  B := simplify(eval(VectorCalculus:-DotProduct(unNorVu, unNorVv)), assume=real);
  C := simplify(eval(VectorCalculus:-DotProduct(unNorVv, unNorVv)), assume=real);
  
  [A, B, C];

end proc:


surTanPln_3D := proc(X, varL, w, s) # X <- Column Vector
  local unXu, unXv; 
  
  unXu := surXuTanVec_3D(X, varL, true);
  unXv := surXvTanVec_3D(X, varL, true);
  
  # Parametric tangent plane to surface X   
  simplify(eval([ X[1] + w*unXu[1] + s*unXv[1], 
                  X[2] + w*unXu[2] + s*unXv[2], 
                  X[3] + w*unXu[3] + s*unXv[3] ]), assume=real);

end proc:  
  
  
surNorPlnXu_3D := proc(X, varL, w, s) # X <- Column Vector
  local unXu, unNorV; 
  
  unXu := surXuTanVec_3D(X, varL, true);
  unNorV := surNorVec_3D(X, varL, true);
  
  # Parametric tangent plane to surface X by vector Xu  
  simplify(eval([ X[1] + w*unXu[1] + s*unNorV[1], 
                  X[2] + w*unXu[2] + s*unNorV[2], 
                  X[3] + w*unXu[3] + s*unNorV[3] ]), assume=real);

end proc:  
  
  
surNorPlnXv_3D := proc(X, varL, w, s) # X <- Column Vector
  local unXv, unNorV; 
  
  unXv := surXvTanVec_3D(X, varL, true);
  unNorV := surNorVec_3D(X, varL, true);
  
  # Parametric tangent plane to surface X by vector Xv   
  simplify(eval([ X[1] + w*unXv[1] + s*unNorV[1], 
                  X[2] + w*unXv[2] + s*unNorV[2], 
                  X[3] + w*unXv[3] + s*unNorV[3] ]), assume=real);

end proc:  
  
  
surGerGeoEqu_3D := proc(X, varL) # X <- Column Vector
  local varU, varV, FFF ,G ,E ,F ,W, 
        Eu, Fu, Gu, Ev, Fv, Gv, eql, eq2;
  
  varU := varL[1];
  varV := varL[2];
  
  FFF  := coeFirFunForm_3D(X, varL);
  
  E  := FFF[1]; 
  F  := FFF[2]; 
  G  := FFF[3];
  
  W  := wFirFunForm_3D(X, varL);
  
  Eu := diff(E, varU); 
  Fu := diff(F, varU); 
  Gu := diff(G, varU); 
  
  Ev := diff(E, varV);
  Fv := diff(F, varV);
  Gv := diff(G, varV);

  eql := simplify(eval( diff(varU(t),t$2) = - subs({varU = varU(t),varV = varV(t)},
                        (G*Eu-F*(2*Fu-Ev))/(2*W))*diff(varU(t),t)^2 
                        - subs({varU = varU(t),varV = varV(t)},(G*Ev-F*Gu)/W)*diff(varU(t),t)* 
                        diff(varV(t),t) - subs({varU = varU(t),varV = varV(t)},
                        (G*(2*Fv-Gu)-F*Gv)/(2*W))*diff(varV(t),t)^2 ), assume=real);
 
  eq2 := simplify(eval( diff(varV(t),t$2) = - subs({varU = varU(t),varV = varV(t)}, 
                        (E*(2*Fu-Ev)-F*Eu)/(2*W))*diff(varU(t),t)^2 
                        - subs({varU = varU(t),varV = varV(t)},(E*Gu-F*Ev)/W)*diff(varU(t),t)*
                        diff(varV(t),t) - subs({varU = varU(t),varV = varV(t)},
                        (E*Gv-F*(2*Fv-Gu))/(2*W))*diff(varV(t),t)^2 ), assume=real);
 
  eql,eq2;

end proc:  
  
  
surTanLine_3D := proc(X, varL, PointL) # X <- Column Vector
  local Xp, up, vp, sTan, sTanp;

  up := PointL[1];
  vp := PointL[2];
  
  Xp := subs(u=up, v=vp, X);
   
  sTan := surTanVec_3D(X, varL, false);
  sTanp := subs(u=up, v=vp, sTan);

  simplify(eval( [Xp[1] + sTanp[1]*t, 
                  Xp[2] + sTanp[2]*t, 
                  Xp[3] + sTanp[3]*t ] ), assume=real);

end proc:


surNorLine_3D := proc(X, varL, PointL) # X <- Column Vector
  local Xp, up, vp, sNor, sNorp;

  up := PointL[1];
  vp := PointL[2];
  
  Xp := subs(u=up, v=vp, X);
   
  sNor := surNorVec_3D(X, varL, false);
  sNorp := subs(u=up, v=vp, sNor);

  simplify(eval( [Xp[1] + sNorp[1]*t, 
                  Xp[2] + sNorp[2]*t, 
                  Xp[3] + sNorp[3]*t ] ), assume=real);

end proc:


Christoffel_3D:= proc(X, varL, indexL) # X <- Column Vector
  local fff, E, F, G, u, v, i, j, k, Eu, Ev, Gu, Gv, Fu, Fv, ChrSym, A; 
  
  fff := coeFirFunForm_3D(X, varL);
     
  E:= fff[1];
  F:= fff[2]; 
  G:= fff[3];
     
  u := varL[1]; #u:= op(1, indets(X,name));
  v := varL[2]; #v:= op(2, indets(X,name));
       
  i:= indexL[1];
  j:= indexL[2];
  k:= indexL[3];
     
  A:= 2*(E*G - F^2);
       
  if i = u and j = u and k = u
  then
    Eu:= diff(E, u);
    Fu:= diff(F, u);
    Ev:= diff(E, v);
				
    ChrSym := (G*Eu - 2*F*Fu + F*Ev)/A
  else
    if i = u and j = u and k = v
    then
      Eu:= diff(E, u);
      Fu:= diff(F, u);
      Ev:= diff(E, v);
				
      ChrSym := (- F*Eu + 2*E*Fu - E*Ev)/A
    else
      if i = u and j = v and k = u
      then
        Ev:= diff(E, v);
        Gu:= diff(G, u);
				
        ChrSym := (G*Ev - F*Gu)/A
      else
        if i = u and j = v and k = v
        then
          Ev:= diff(E, v);
          Gu:= diff(G, u);
				    
          ChrSym := (E*Gu - F*Ev)/A
        else
          if i = v and j = v and k = u
          then
            Fv:= diff(F, v);
            Gu:= diff(G, u);	
            Gv:= diff(G, v);
						  
            ChrSym := (- F*Gv + 2*G*Fv - G*Gu)/A
          else
            if i = v and j = v and k = v
            then
              Fv:= diff(F, v);
              Gu:= diff(G, u);	
              Gv:= diff(G, v);
				    
              ChrSym := (E*Gv - 2*F*Fv + F*Gu)/A
            end if;
          end if;
        end if;
      end if;
    end if;
  end if; 

  simplify(eval( ChrSym ), assume=real);
     
end proc:



















  
